# 加载必要的库
library(sf)
library(tmap)
library(gridExtra)
library(purrr)
library(grid)
library(dplyr)


source("./code/2.3 人口数估发病率.R")

# 逐年_空间自相关------

# 假设 combined_df 和 map_data 已经加载到环境中
combined_df <- readxl::read_excel( "./output/表格3.4.2_LISA_历年.xlsx") 
is.na(combined_df$Pr_Ii) %>% sum()
# combined_df %>% View()

# combined_df <- combined_df %>%  #20250325 不需要设定是否为缺失；
#   mutate(Pr_Ii = ifelse(is.na(Pr_Ii),
#                         0, Pr_Ii))
# nameCE_county$county_code  <-  as.character(nameCE_county$county_code)
map_data <- st_read("./data/origin/新疆维吾尔自治区.json") %>%
  st_transform(crs = 4326) %>%
  mutate(NAME = name, GB1999 = as.character(adcode), city_code = substr(GB1999,1,4)) 

# 检查哪些几何对象是无效的
invalid_index <- which(!st_is_valid(st_geometry(map_data)))

if (length(invalid_index) > 0) {
  cat("发现无效几何对象，尝试修复...\n")
  
  # 对无效几何对象进行修复
  map_data$geometry[invalid_index] <- st_make_valid(st_geometry(map_data)[invalid_index])
  
  # 再次检查是否仍然有无效几何对象
  invalid_after_first_fix <- which(!st_is_valid(st_geometry(map_data)))
  
  if (length(invalid_after_first_fix) > 0) {
    cat("仍有无效几何对象，尝试使用缓冲处理修复...\n")
    
    # 使用 st_buffer 修复剩余的无效几何对象
    buffer_dist <- 0.000001  # 使用一个非常小的正值避免重复边问题
    map_data$geometry[invalid_after_first_fix] <- st_buffer(
      st_geometry(map_data)[invalid_after_first_fix],
      dist = buffer_dist
    )
    
    # 再次检查是否仍然有无效几何对象
    invalid_after_buffer <- which(!st_is_valid(st_geometry(map_data)))
    
    if (length(invalid_after_buffer) > 0) {
      cat("仍有无效几何对象，尝试简化几何对象...\n")
      
      # 尝试简化几何对象
      simplify_tolerance <- 0.0001  # 根据需要调整简化容差
      map_data$geometry[invalid_after_buffer] <- st_simplify(
        st_geometry(map_data)[invalid_after_buffer],
        dTolerance = simplify_tolerance
      )
      
      # 最后再检查一次有效性
      final_invalid <- which(!st_is_valid(st_geometry(map_data)))
      if (length(final_invalid) > 0) {
        warning("仍有无效几何对象无法自动修复，请手动检查这些对象。\n")
      } else {
        cat("所有几何对象已成功修复。\n")
      }
    } else {
      cat("所有几何对象已通过缓冲处理成功修复。\n")
    }
  } else {
    cat("所有几何对象已通过 st_make_valid 成功修复。\n")
  }
} else {
  cat("所有几何对象都是有效的，无需修复。\n")
}

# 确保 map_data 的几何列是最新的
map_data <- st_as_sf(map_data)

# 检查哪些几何对象是无效的
invalid_index <- which(!st_is_valid(st_geometry(map_data)))

# 打印结果以确认
print(head(map_data))

# 合并空间数据和属性数据
map_data <- merge(map_data, combined_df, by.x = "adcode", by.y = "county_code")

names(map_data)


# 创建一个列表来存储逐年地图
yearly_maps <- list()

# 循环绘制每年的空间自相关地图（2005-2023）
for (year in 2005:2023) {
  yearly_data <- map_data[map_data$year == year, ] %>% 
    mutate(cluster_type = case_when(
      Pr_Ii >= 0.05 ~ "Not Significant",
      TRUE ~ cluster_type
    )) 
  
  # 定义颜色映射
  cluster_colors <- c(
    "High-High" =  "#d73027", 
    "High-Low" = "#fc8d59",
    "Low-High" = "#91bfdb", 
    "Low-Low" =  "blue",#"#4575b4",
    "Not Significant" = "white"# Other category
  )
  
  # 绘制地图
  yearly_map <- tm_shape(yearly_data) +
    tm_polygons(
      col = "cluster_type",
      palette = cluster_colors,
      title = "Cluster Type",
      style = "cat",
      border.col = "gray40",
      border.lwd = 0.4,
      alpha = 0.9,
      legend.show = F
    ) +
    tm_layout(
      main.title = paste0(year),
      # main.title.size = 1.2,
      main.title.position = "center",
      main.title.fontface = "bold",
      fontfamily = "Arial",
      legend.position = c("right", "bottom"),
      legend.bg.color = "white",
      legend.bg.alpha = 0.8,
      legend.frame = FALSE,
      # legend.text.size = 0.85,
      # legend.title.size = 1.0,
      legend.show = F,
      frame = FALSE,
      outer.margins = c(0.02, 0, 0.02, 0)
    ) #+
  tm_compass(position = c("left", "top"), size = 2) +  # 添加指北针
    tm_scale_bar()
  
  # 将每年的地图存入列表
  yearly_maps[[as.character(year)]] <- yearly_map
}

yearly_maps[[2]]

# 定义颜色映射
cluster_colors <- c(
  "High-High" =  "#d73027", 
  "High-Low" = "#fc8d59",
  "Low-High" = "#91bfdb", 
  "Low-Low" =  "blue",#"#4575b4",
  "Not Significant" = "white"
)

# 创建第20张图：比例尺、指北针和图例

legend_map <- tm_shape(map_data) +
  tm_polygons() +
  tm_layout(bg.color = "transparent") +
  tm_compass(#position = c("left", "top"),
    position = c(0.2,0.8),
    size = 1.5) +
  tm_scale_bar(
    # position = c("right", "top"),
    position = c(0.2,0.65),
    breaks = c(0,400,800),
    text.size = 0.7) +
  tm_add_legend(
    type = "fill",
    labels = c("High-High",  "High-Low", "Low-High","Low-Low", "Not Significant"),
    col = c("#d73027", "#fc8d59",  "#91bfdb","blue","white"),
    title = "Cluster Category",
    border.lwd = 0.5
  ) +
  tm_layout(
    legend.only = TRUE,
    legend.position = c(0.2, 0.2),
    legend.width = 1.5,
    legend.title.size = 1.2,
    legend.text.size = 0.8
  )
legend_map

# 将 yearly_maps 和 legend_map 合并为一个列表
combined_list <- c(yearly_maps, list(legend_map))

# 组合输出 (5列4行布局)
final_map <- tmap_arrange(
  combined_list,
  ncol = 5,  # 每行 5 张图
  nrow = 4,  # 每列 4 张图
  outer.margins = 0.03
  # # 关键参数优化
  # outer.margins = c(0, 0, 0, 0),  # 上下左右外边距归零
  # asp = NULL,                     # 禁用固定长宽比
  # widths = rep(1, 5),             # 列宽均等分配
  # heights = rep(1, 4),            # 行高均等分配
  # padding = unit(0.1, "cm")       # 子图间距调小
  # asp = 0.6 # 要注释掉，不然图很小、图间空白很大
)

final_map

tmap_save(
  final_map,
  filename = "./主图/LISA地图/逐年_LISA图.png",
  width = 45,     # 宽度 45cm（期刊通常接受 10-50cm）
  height = 35,    # 高度 35cm
  units = "cm",   # 单位厘米（期刊推荐）
  dpi = 600       # 分辨率 600 DPI（期刊要求 ≥ 300 DPI）
)




# 
# # 逐年月_空间自相关------
# 
# # 过滤出2023年的数据
# rate_month
# 
# year_2023_data <- yearly_incidence_county %>%
#   filter(year == 2023) %>%
#   select(county_code, incidence_rate)
# 
# # 确保所有区县都在数据中，并用0填充缺失值
# year_2023_map_data <- map_data %>%
#   left_join(year_2023_data, by = c("GB1999" = "county_code")) %>%
#   mutate(incidence_rate = replace_na(incidence_rate, 0))
# 
# # 创建邻居关系
# nb_q <- poly2nb(st_geometry(map_data), queen = TRUE)
# 
# # 创建权重列表，并明确指定 zero.policy = TRUE
# lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)
# length(lw_q$neighbours)
# 
# # 执行全局莫兰检验
# global_moran_result <- moran.test(year_2023_map_data$incidence_rate, listw = lw_q,
#                                   randomisation = FALSE, zero.policy = TRUE)
# 
# # 查看全局莫兰检验结果
# print(global_moran_result)
# 
# # LISA结果
# local_moran_result <- localmoran(year_2023_map_data$incidence_rate, listw = lw_q,
#                                  zero.policy = TRUE)
# 
# # 整理局部莫兰检验结果为数据框
# local_moran_df <- data.frame(
#   county_code = year_2023_map_data$GB1999,
#   name = year_2023_map_data$name,
#   incidence_rate = year_2023_map_data$incidence_rate,
#   Ii = local_moran_result[, 1],
#   E_Ii = local_moran_result[, 2],
#   Var_Ii = local_moran_result[, 3],
#   Z_Ii = local_moran_result[, 4],
#   Pr_Ii = local_moran_result[, 5]
# )
# 
# # 查看局部莫兰检验结果的前几行
# print(head(local_moran_df))
# 
# 
# 
# 
# ## 通义千问代码模式----
# library(dplyr)
# library(spdep)
# library(sf)
# 
# # 假设map_data已经加载并且包含地理信息，并且GB1999字段对应county_code
# 
# # 创建邻居关系
# nb_q <- poly2nb(st_geometry(map_data), queen = TRUE)
# 
# # 创建权重列表，并明确指定 zero.policy = TRUE
# lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)
# 
# # 定义一个函数来计算每年每月的LISA
# calculate_monthly_lisa <- function(data) {
#   # 确保所有区县都在数据中，并用0填充缺失值
#   map_data_with_rate <- map_data %>%
#     left_join(data, by = c("GB1999" = "地区编码")) %>%
#     mutate(rate = replace_na(rate, 0))
#   
#   # 执行全局莫兰检验
#   global_moran_result <- moran.test(map_data_with_rate$rate, listw = lw_q,
#                                     randomisation = FALSE, zero.policy = TRUE)
#   
#   # LISA结果
#   local_moran_result <- localmoran(map_data_with_rate$rate, listw = lw_q,
#                                    zero.policy = TRUE)
#   
#   # 整理局部莫兰检验结果为数据框
#   local_moran_df <- data.frame(
#     地区编码 = map_data_with_rate$GB1999,
#     name = map_data_with_rate$name,
#     rate = map_data_with_rate$rate,
#     Ii = local_moran_result[, 1],
#     E_Ii = local_moran_result[, 2],
#     Var_Ii = local_moran_result[, 3],
#     Z_Ii = local_moran_result[, 4],
#     Pr_Ii = local_moran_result[, 5],
#     year = data$year[1],
#     month = data$month[1]
#   )
#   
#   return(local_moran_df)
# }
# 
# # 对每个年份和月份分组并应用上述函数
# monthly_lisa_results <- rate_month %>%
#   group_by(year, month) %>%
#   do(calculate_monthly_lisa(.)) %>%
#   ungroup()
# # 
# # # 查看结果
# # print(head(monthly_lisa_results))
# # monthly_lisa_results %>% dim()
# # monthly_lisa_results$year %>% unique()
# # 
# # writexl::write_xlsx(monthly_lisa_results, "./output/表格3.4.2_LISA_历年月.xlsx")
# # 
# # 
# # # 初始化结果列表
# # cluster_summary <- list()
# # 
# # # 整理局部莫兰检验结果
# # for (yr_mon in unique(paste(monthly_lisa_results$year, monthly_lisa_results$month))) {
# #   yearly_monthly_local_moran <- monthly_lisa_results %>%
# #     filter(paste(year, month) == yr_mon) %>%
# #     filter(Pr_Ii < 0.05)
# #   
# #   yearly_monthly_local_moran <- yearly_monthly_local_moran %>%
# #     mutate(cluster_type = case_when(
# #       Z_Ii > 0 & Ii > E_Ii ~ "High-High",
# #       Z_Ii < 0 & Ii < E_Ii ~ "Low-Low",
# #       Z_Ii > 0 & Ii < E_Ii ~ "High-Low",
# #       Z_Ii < 0 & Ii > E_Ii ~ "Low-High"
# #     ))
# #   
# #   cluster_counts <- yearly_monthly_local_moran %>%
# #     group_by(cluster_type) %>%
# #     summarise(count = n(), 
# #               counties = paste(name, collapse = ", ")) %>%
# #     ungroup() %>%
# #     complete(cluster_type = c("High-High", "Low-Low", "High-Low", "Low-High"), fill = list(count = 0, counties = ""))
# #   
# #   cluster_summary[[yr_mon]] <- cluster_counts
# # }
# # 
# # 
# # # 查看结果
# # print(str(cluster_summary))
# 
# 
# # 初始化一个空列表来保存调整后的数据框
# adjusted_dataframes <- list()
# 
# # 遍历 monthly_lisa_results 数据框，为每个数据框添加 'cluster_type' 列，并存入 adjusted_dataframes 列表
# for (yr_mon in unique(paste(monthly_lisa_results$year, monthly_lisa_results$month))) {
#   df <- monthly_lisa_results %>%
#     filter(paste(year, month) == yr_mon) %>%
#     filter(Pr_Ii < 0.05)
#   
#   df <- df %>%
#     mutate(cluster_type = case_when(
#       Z_Ii > 0 & Ii > E_Ii ~ "High-High",
#       Z_Ii < 0 & Ii < E_Ii ~ "Low-Low",
#       Z_Ii > 0 & Ii < E_Ii ~ "High-Low",
#       Z_Ii < 0 & Ii > E_Ii ~ "Low-High"
#     ))
#   
#   adjusted_dataframes[[yr_mon]] <- df
# }
# 
# # 使用 bind_rows() 将所有调整后的数据框合并成一个大的数据框
# combined_df <- bind_rows(adjusted_dataframes)
# 
# # 查看结果
# head(combined_df)
# 
# # 保存结果到Excel文件
# writexl::write_xlsx(combined_df, "./output/表格3.4.2_LISA_历年月.xlsx")
# 
# # 整理局部莫兰检验结果
# cluster_summary <- list()
# 
# # 整理局部莫兰检验结果
# for (yr_mon in unique(paste(combined_df$year, combined_df$month))) {
#   yearly_monthly_local_moran <- combined_df %>%
#     filter(paste(year, month) == yr_mon) %>%
#     filter(Pr_Ii < 0.05)
#   
#   cluster_counts <- yearly_monthly_local_moran %>%
#     group_by(cluster_type) %>%
#     summarise(count = n(), 
#               counties = paste(name, collapse = ", ")) %>%
#     ungroup() %>%
#     complete(cluster_type = c("High-High", "Low-Low", "High-Low", "Low-High"), fill = list(count = 0, counties = ""))
#   
#   cluster_summary[[yr_mon]] <- cluster_counts
# }
# 
# # 合并结果
# cluster_summary_df <- bind_rows(cluster_summary, .id = "year_month")
# 
# # 调整列顺序以便查看
# cluster_summary_df <- cluster_summary_df %>%
#   separate(year_month, into = c("year", "month"), sep = " ") %>%
#   select(year, month, cluster_type, count, counties)
# 
# # 查看局部莫兰检验结果
# print(head(cluster_summary_df))
# 
# # 步骤 5：保存结果
# writexl::write_xlsx(cluster_summary_df, "./output/表格3.4.2_局部莫兰_LISA_计数_逐年月.xlsx")
# 
# # 步骤 6：将局部莫兰指数转宽数据
# cluster_summary_df1 <- cluster_summary_df %>% dplyr::select(-counties)
# cluster_summary_df2 <- cluster_summary_df1 %>% 
#   pivot_wider(names_from = cluster_type, values_from = count, values_fill = 0) 
# 
# # 查看宽数据
# print(head(cluster_summary_df2))
# 
# # 保存宽数据到Excel文件
# writexl::write_xlsx(cluster_summary_df2, "./output/表格3.4.3_LISA_宽数据_历年月.xlsx")
# 
# # 步骤 7：将历年聚集类型绘图
# 
# cluster_summary_df2$year <- as.numeric(cluster_summary_df2$year)
# cluster_summary_df2$month <- as.numeric(cluster_summary_df2$month)
# 
# # 创建日期变量
# cluster_summary_df2$date <- as.Date(paste(cluster_summary_df2$year, cluster_summary_df2$month, "01", sep = "-"))
# 
# # 折线图：不同类型聚集区县数量
# ggplot(cluster_summary_df2, aes(x = date)) +
#   geom_line(aes(y = `High-High`, color = "High-High"), size = 1) + 
#   geom_point(aes(y = `High-High`, color = "High-High"), size = 3) + 
#   geom_line(aes(y = `Low-Low`, color = "Low-Low"), size = 1) + 
#   geom_point(aes(y = `Low-Low`, color = "Low-Low"), size = 3) + 
#   geom_line(aes(y = `High-Low`, color = "High-Low"), size = 1) + 
#   geom_point(aes(y = `High-Low`, color = "High-Low"), size = 3) + 
#   geom_line(aes(y = `Low-High`, color = "Low-High"), size = 1) + 
#   geom_point(aes(y = `Low-High`, color = "Low-High"), size = 3) + 
#   labs(
#     title = "Local Moran's I Analysis of Brucellosis Incidence in Xinjiang (2005-2023)",
#     x = "Date",
#     y = "Number of Counties",
#     color = "Cluster Type"
#   ) +
#   theme_minimal() + # 使用简洁的主题
#   scale_color_manual(values = RColorBrewer::brewer.pal(4, "Set3")) + # 使用 Set3 色系
#   theme(
#     plot.title = element_text(hjust = 0.5), # 居中标题
#     axis.text.x = element_text(angle = 45, hjust = 1) # 调整x轴标签角度以便阅读
#   )
# 
# ggsave("./figs/图3.4.2 折线图_LISA_历年.png")
# 
# cluster_summary_df_month <- cluster_summary_df %>% 
#   group_by(year, cluster_type) %>% summarise(count = sum(count))
# 
# 
# cluster_summary_df_month$year <- as.numeric(cluster_summary_df_month$year)
# 
# # 折线图：不同类型聚集区县数量
# ggplot(cluster_summary_df_month, aes(x = year, y = count, 
#                                      color = cluster_type, group = cluster_type)) +
#   geom_line(size = 1) + # 添加线条
#   geom_point(size = 3) + # 添加点，使每个数据点更明显
#   labs(
#     title = "Local Moran's I Analysis of Brucellosis Incidence in Xinjiang (2005-2023)",
#     x = "Year",
#     y = "Number of Counties",
#     color = "Cluster Type"
#   ) +
#   theme_minimal() + # 使用简洁的主题
#   scale_color_manual(values = RColorBrewer::brewer.pal(4, "Set3")) + # 使用 Set3 色系
#   theme(
#     plot.title = element_text(hjust = 0.5), # 居中标题
#     axis.text.x = element_text(angle = 45, hjust = 1) # 调整x轴标签角度以便阅读
#   )
# 
# ggsave("./figs/图3.4.3 折线图_LISA_历年月.png")
